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1. INTRODUCTION 

We have earlier in a sequence of papers [Logg 2003a; 2003b; 2006] introduced 
the multi-adaptive Galerkin methods mcG(g) and mdG(g) for the approximate 
(numerical) solution of ODEs of the form 

u(*) = /(«(t),t), 1 6(0,1*], 
u(0) = n , 
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where u : [0,T] — > R N is the solution to be computed, u € R N a given initial 
value, T > a given final time, and / : R N x (0, T] — > R N a given function that is 
Lipschitz continuous in u and bounded. 

The multi-adaptive Galerkin methods mcG(g) and mdG(g) extend the standard 
mono- adaptive continuous and discontinuous Galerkin methods cG{q) and dG(g), 
studied before in [Hulme 1972b; 1972a; Jamet 1978; Delfour et al. 1981; Eriksson 
et al. 1985; Johnson 1988; Eriksson and Johnson 1991; 1995a; 1995b; 1995c; Eriks- 
son et al. 1998; Eriksson et al. 1995; Estep 1995; Estep and French 1994; Estep et al. 
2000; Estep and Williams 1996; Estep and Stuart 2002], by allowing individual time 
step sequences h — ki(t) for the different components Ui = Ui(t), i = 1,2, ... ,N, 
of the approximate solution U ~ u of the initial value problem (1). For related 
work on local time-stepping, see also [Hughes et al. 1983a; 1983b; Makino and 
Aarseth 1992; Dave et al. 1997; Alexander and Agnor 1998; Osher and Sanders 
1983; Flaherty et al. 1997; Dawson and Kirby 2001; Lew et al. 2003; Engstlcr and 
Lubich 1997; Savcenco et al. 2005; Savcenco 2008]. In comparison with existing 
method for local time-stepping, the main advantage of the multi-adaptive Galerkin 
methods mcG(g) and mdG(g) is the automatic local step size selection based on a 
global a posteriori error estimate built into these methods. 

In the current paper, we discuss important aspects of the implementation of 
multi-adaptive Galerkin methods. While earlier results on multi-adaptive time- 
stepping presented in [Logg 2003a; 2003b; 2006] include the formulation of the 
methods, a priori and a posteriori error estimates, together with a proof-of-concept 
implementation and results for a number of model problems, the current paper ad- 
dresses the important issue of efficiently implementing the multi-adaptive methods 
with minimal overhead as compared to standard mono-adaptive solvers. For many 
problems, in particular when the propagation of the solution is local in space and 
time, the potential speedup of multi-adaptivity is large, but the actual speedup may 
be far from the ideal speedup if the overhead of the more complex implementation 
is significant. 

1.1 Implementation 

The algorithms presented in this paper are implemented by the multi-adaptive 
ODE-solver available in DOLFIN [Logg et al. ; Hoffman and Logg 2002], Dynamic 
Object-oriented Library for FINitc clement computation, which is the C++ inter- 
face of the new open-source software project FEniCS [FEniCS 2008; Logg 2007; 
Dupont et al. 2003] for the automation of Computational Mathematical Modeling 
(CMM). The multi- adaptive solver in DOLFIN is based on the original implemen- 
tation Tanganyika, presented in [Logg 2003b], but has been completely rewritten 
for DOLFIN and is actively developed by the authors. 

1.2 Obtaining the software 

DOLFIN is licensed under the GNU (Lesser) General Public License [Free Software 
Foundation 1999], which means that anyone is free to use or modify the software, 
provided these rights are preserved. The complete source code of DOLFIN, includ- 
ing numerous example programs, is available at the DOLFIN web page [Logg et al. 
]• 
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1.3 Notation 

The following notation is used throughout this paper: Each component Ui(t), 
i = 1,...,N, of the approximate m(c/d)G(g) solution U(t) of (1) is a piecewise 
polynomial on a partition of (0, T] into m, sub-intervals. Sub-interval j for compo- 
nent i is denoted by hj = (tij-i,tij], and the length of the sub-interval is given by 
the local time step kij — t^ —ti,j-\. We shall sometimes refer to Iij as an element. 
This is illustrated in Figure 1. On each sub- interval Iij, Ui\i tj is a polynomial of 
degree at most qij. 

Furthermore, we shall assume that the interval (0,T] is partitioned into blocks 
between certain synchronized time levels = T < Ti < . . . < T M = T. For 
each T n , n = 0,1,..., M and each i = 1,2, ... ,N, we require that there is a 
< j < mi such that tij = T n . We refer to the collection of local intervals between 
two synchronized time levels T„_i and T„ as a time slab. We denote the length of 
a time slab by K n = T n — T n _i. 



t 

T 



Tn— 1 Kn T n 



Fig. 1. Individual partitions of the interval (0,T] for different components. Elements between 
common synchronized time levels are organized in time slabs. In this example, we have N = 6 
and M = 4. 



1.4 Outline of the paper 

We first give an introduction to multi-adaptive time-stepping in Section 2. We then 
present the key algorithms used by the multi-adaptive ODE solver of DOLFIN in 
Section 3, followed by a discussion of data structures for efficient representation and 
interpolation of multi-adaptive solutions in Section 4. In Section 5, we discuss the 
efficiency of multi-adaptive time-stepping and in Section 6, we present a number 
of numerical examples that demonstrate the efficiency of the proposed algorithms 
and data structures. Finally, we give some concluding remarks in Section 7. 
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2. MULTI-ADAPTIVE TIME-STEPPING 

In this section, we give a quick introduction to multi-adaptive time-stepping, in- 
cluding the formulation of the methods, error estimates and adaptivity. For a more 
detailed account, we refer the reader to [Logg 2003a; 2003b; 2006]. 

2.1 Formulation of the methods 

The mcG(g) and mdG(q) methods are obtained by multiplying the system of equa- 
tions (1) with a suitable test function v, to obtain the following variational problem: 
Find U e V with 17(0) = u , such that 

(v,U)dt = [ (v,f(U,-))dt VveV, (2) 
Jo 

where (•, •) denotes the standard inner product on R w and (V, V) is a suitable pair 
of discrete function spaces, the test and trial spaces respectively. 

For the standard cG(g) method, the trial space V consists of the space of con- 
tinuous piecewise polynomial vector-valued functions of degree q = q(t) on a par- 
tition = to < t\ < ••• < t M = T and the test space V consist of the space 
of (possibly discontinuous) piecewise polynomial vector-valued functions of degree 
q — 1 on the same partition. The multi-adaptive mcG(g) method extends the stan- 
dard cG(q) method by extending the test and trial spaces to piecewise polynomial 
spaces on individual partitions of the time interval that satisfy the constraints in- 
troduced in the previous section and illustrated in Figure 1. Thus, each component 
Ui = Ui(t) is continuous and a piecewise polynomial on the individual partition 
= Uo < tn < ■ ■ ■ < t imi = T for i = 1, 2, . . . , N. 

For the standard dG(g) method, the test and trial spaces are equal and consist of 
the space of (possibly discontinuous) piecewise polynomial vector- valued functions 
of degree q = q{t) on a partition = to < t\ < • • • < tia = T, which extends 
naturally to the multi- adaptive mdG(g) method by allowing each component of 
the test and trial functions to be a piecewise polynomial on its own partition of 
the time interval as above for the mcG(q) method. Note that for both the dG(q) 
method and the mdG(q) method, the integral / T (v, U) dt in (2) must be treated 
appropriately at the points of discontinuity, see [Logg 2003a] . 

Both in the case of the mcG(g) and mdG(g) methods, the variational problem (2) 
gives rise to a system of discrete equations by expanding the solution U in a suitable 
basis on each local interval , 

Hi 

= y ' (,ijm<ftijm, (3) 

where {£,ij m }m=o are the degrees of freedom for Ui on 1^ and {<frijm}m=o 1S a 
suitable basis for P qij (7^ ) . For any particular choice of quadrature, the resulting 
system of discrete equations takes the form of an implicit Runge-Kutta method on 
each local interval . The discrete equations take the form 

Hijrn = e« + % E «&« ] W{T^(S^)),T^{S^)), (4) 
n=0 
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for m = 0, . . . , qij, where {C }m=o, n =o are weights, T y maps to (0, 1], r y (t) = 

o are quadrature points defined on [0,1]. 
Note that we have here assumed that the number of quadrature points is equal to 
the number of nodal points. See [Logg 2003a] for a discussion of suitable quadrature 
rules and basis functions. 

2.2 Error estimates and adaptivity 

The global error e = U — u of the approximate solution U of (1) may be bounded 
in terms of computable quantities. Such an a posteriori error estimate is proved 
in [Logg 2003a], both for the mcG(g) and mdG(g) methods. The a posteriori error 
estimate provides a bound for any given linear functional M. : R N — > M of the 
global error e(T) at the final time, such as the error e^(T) in a single component. 
Bounds for the error itself in various norms may also be approximated. Below, we 
state the basic a posteriori error estimate for the mcG(g) method and refer to [Logg 
2003a] for a complete discussion, including error estimates for mdG(g). 
For the mcG(<;) method, the error estimate takes the following form: 

N 

\M(e(T))\<E = Y, S i( T )™%{ C t k T\Ri\}, (5) 

i—l 

Here, R = U — f(U, •) denotes the residual of the computed solution, C\ = Ci(t) 
denotes an interpolation constant (which may be different for each local interval) 
and Si (T) denotes a stability factor that measures the rate of propagation of local 
errors for component Ui (the influence of a nonzero residual in component Ui on the 
size of the error in the given functional) . By selecting the local time steps ki = ki (t) 
such that E = TOL for a given tolerance TOL, one may thus guarantee that the 
error in the functional Ai is bounded by the given tolerance, |jM(e(T))| < TOL. 

Comparing to standard Runge-Kutta methods for the solution of initial value 
problems, the stability factor quantifies the relationship between the "local error" 
and the global error. Note that alternatively, the stability information may be 
kept as a local time-dependent stability weight for more fine-grained control of the 
contributions to the global error. The stability factors are obtained by solving a dual 
problem of (1) for the given functional A4, see [Eriksson et al. 1995; Logg 2003a]. 
The particular form of the dual problem for (1) will be discussed in Section 3.5. 

The individual time steps may be chosen so as to equidistribute the error in the 
different components in an attempt to satisfy 

Cijk'% max \Ri\ - TOL/(JVSj(T)), (6) 

for each local time interval Iij. This may be done in an iterative fashion, as outlined 
in the following basic adaptive algorithm: 

(0) Assume S l (T) = 1 for i = 1, 2, . . . , N; 

(1) Solve the primal problem with time steps based on (6); 

(ii) Solve the dual problem and compute the stability factors; 

(iii) Compute an error bound E based on (5); 

(iv) If E < TOL then stop; if not go back to (i). 
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3. ALGORITHMS 

We present below a collection of key algorithms for multi-adaptive time-stepping. 
The algorithms are given in pseudo-code and where appropriate we give remarks 
on how the algorithms have been implemented in C++ for DOLFIN. In most cases, 
we present simplified versions of the algorithms with focus on the most essential 
steps. 

3.1 General algorithm 

The general multi-adaptive time-stepping algorithm is Algorithm 1. Starting at 
t = 0, the algorithm creates a sequence of time slabs until the given end time T 
is reached. In each macro time step, Algorithm 2 (CreateTimeSlab) is called to 
create a time slab covering an interval [T n _ 1 ,T n ] such that T n < T. For each time 
slab, the system of discrete equations is solved iteratively, using direct fixed-point 
iteration or a preconditioned Newton's method, until the discrete equations given 
by the mcG(g) or mdG(g) method have converged. 



Algorithm 1 U — Integratc(ODE) 
t <- 

while t < T 

{time slab, t} <- CreateTimeSlab({l, ...,N},t, T) 
SolvcTimeSlab(time slab) 
end while 



The basic forward integrator, Algorithm 1, can be used as the main component 
of an adaptive algorithm with automated error control of the computed solution 
as outlined in Section 2. In each iteration, the primal problem (1) is solved using 
Algorithm 1. An ODE of the form (1) representing the dual problem is then created 
and solved using Algorithm 1. It is important to note that both the primal and 
the dual problems may be solved using the same algorithm, but with (possibly) 
different time steps, tolerances, methods, and orders. When the solution of the 
dual problem has been computed, the stability factors {Si{T)}f =1 and the error 
estimate may be computed. 

3.2 Recursive construction of time slabs 

In each step of Algorithm 1, a new time slab is created between two synchronized 
time levels T„_i and T n . The time slab is organized recursively as follows. The 
root time slab covering the interval [T n _i, T n ] contains a non-empty list of elements, 
which we refer to as an element group, and a possibly empty list of time slabs, which 
in turn may contain nested groups of elements and time slabs. Each such element 
group together with the corresponding nested set of element groups is referred to 
as a sub-slab. This is illustrated in Figure 2. 

To create a time slab, we first compute the desired time steps for all components 
as given by the a posteriori error estimate (5). We discuss in detail the time step 
selection below in Section 3.4. A threshold 9K is then computed based on the 
maximum time step K among the components and a fixed parameter 6 <E (0, 1) 
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Fig. 2. The recursive organization of the time slab. Each time slab contains an element group 
and a list of recursively nested time slabs. The root time slab in the figure contains one element 
group of one element and three sub-slabs. The first of these sub-slabs contains an element group 
of two elements and two nested sub-slabs, and so on. The root time slab recursively contains a 
total of nine element groups and 33 elements. 



controlling the density of the time slab. The components are partitioned into two 
sets based on the threshold, and a large time step K_ is selected to be the smallest 
time step among the components in the set with large time steps as described in 
Algorithm 3 and illustrated in Figure 3. For each component in the group with 
large time steps, an element is created and added to the element group of the time 
slab. The remaining components with small time steps are processed by a recursive 
application of this algorithm for the construction of time slabs. 

We organize the recursive construction of time slabs as described by Algorithms 
2, 3, 4, and 5. The recursive construction simplifies the implementation; each 
recursively nested sub-slab can be considered as a sub-system of the ODE. Note 
that the element group containing elements for components in group I\ is created 
before the recursively nested sub-slabs for components in group in- The tree of 
time slabs is thus created recursively breadth- first, which means in particular that 
the element for the component with the largest time step is created first. 

Algorithm 3 for the partition of components can be implemented efficiently using 
the function std: : partition () , which is part of the Standard CH — h Library. 

3.3 Solving the system of discrete equations 

On each time slab T n , n = 1,2,..., M, we need to solve a system of equations for 
the degrees of freedom on the time slab. On each local interval Uj 6 T n , these 
equations are given by (4). Depending on the properties of the given system (1), 
different solution strategies for the time slab system (4) may be appropriate as 
outlined below. 
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Algorithm 2 {time slab, T n } = CrcateTimcSlab(componcnts, T„_i, T) 

{Io, I\, K} Partition(components) 
if T n ^+K <T 
T„ <- T„_i + K 

else 
end if 

element group CreateElements(ii, T„_i, T„) 
time slabs <— CreateTimeSlabs(io, T n _\, T n ) 
time slab {clement group, time slabs} 



Algorithm 3 {/o, Ii, K} — Partition(components) 
Io ^0 

Ji ^0 

if <— maximum time step within components 
for each component 

k <— time step of component 

if k < 6K 

la h U {component} 

else 

/i<-/iU {component} 
endif 
end for 

K_ minimum time step within I\ 



Algorithm 4 elements = CrcatcElements(components, T„_i, T„) 

elements «— 

for each component 

create element for component on [T n -i,T n ] 

elements <— elements U element 
end for 



Algorithm 5 time slabs = CreateTimeSlabs(components, T„_i, T n ) 

time slabs <— 
* <- T„_i 
while f < T 

{time slab, £} «— CreateTimcSlab(components, i, T„) 
time slabs ■(— time slabs U time slab 
end while 
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BK K K 
Fig. 3. The partition of components into groups of small and large time steps for 9 = 1/2. 



3.3.1 Direct fixed-point iteration. In the simplest case, the time slab system is 
solved by direct fixed-point iteration on (4) for each element in the time slab. The 
fixed-point iteration is performed in a forward fashion, sweeping over the elements 
in the time slab in the same order as they are created by Algorithm 2. In particular, 
this means that for each component in the time slab system, the end-time value on 
each element is updated before the degrees of freedom for the following element. 
Thus, for each element 6 T n , we compute the degrees of freedom {£ym}f=o 
according to 

la 

tiim = ^0 + h J2 W ™ fi(U^\s^)),Tv\8^)\ TU = 0, 1, . . . , Qij. (7) 

71=0 

Direct fixed-point iteration converges if the system is non-stiff and typically only a 
few iterations are needed. In fact, one may consider a system to be stiff if direct 
fixed-point iteration does not converge. 

3.3.2 Damped fixed-point iteration. If the system is stiff, that is, direct fixed- 
point iteration does not converge, one may introduce a suitable amount of damping 
to adaptively stabilize the fixed-point iteration. The fixed-point iteration (7) may 
be written in the form 

where £ is the vector of degrees of freedom for the solution on the time slab. We 
modify the fixed-point iteration by introducing a damping parameter a: 
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In [Logg 2004], a number of different strategies for the selection of the damping 
parameter a are discussed. We mention two of these strategies here. The first 
strategy chooses a based on the diagonal derivatives dfi/dui, i = 1, 2, . . . , N, cor- 
responding to a modified Newton's method where the Jacobian is approximated by 
a diagonal matrix. This strategy works well for systems with a diagonally domi- 
nant Jacobian, including many systems arising when modeling chemical reactions. 
The second strategy adaptively chooses a scalar a based on the convergence of the 
fixed-point iterations. 

3.3.3 Newton's method. Alternatively, one may apply Newton's method directly 
to the full system of equations (7) associated with each time slab. The linear system 
in each Newton iteration may then be solved either by a direct method or an 
iterative method such as a Krylov subspace method in combination with a suitable 
preconditioner, depending on the characteristics of the underlying system (1). In 
addition, one may also apply a special preconditioner that improves the convergence 
by propagating values forward in time within the time slab. Note that if the multi- 
adaptive efficiency index is large (see Section 5 below), then the time slab system is 
not significantly larger than the corresponding time slab system for a mono-adaptive 
method. 

3.3.4 Choosing a solution strategy. Ultimately, an intelligent solver should au- 
tomatically choose a suitable algorithm for the solution of the time slab system. 
Thus, the solver may initially try direct fixed-point iteration. If the system is stiff, 
the solver switches to adaptive fixed-point iteration (as outlined in [Logg 2004]). 
Finally, if the adaptive fixed-point iteration converges slowly, the solver may switch 
to Newton's method. 

3.3.5 Interpolation of the solution. To update the degrees of freedom on an 
element according to (7), the appropriate component fi of the right-hand side of 
(1) needs to be evaluated at the set of quadrature points. In order for /, to be 
evaluated, each component \Jy of the computed solution U on which /, depends, 
needs to be evaluated at the quadrature points. We let Si C {1, ...,JV} denote 
the sparsity pattern of component t/j, that is, the set of components on which fi 
depends, 

5, = {*' e {1, • . • , N) : dfi/dm, + 0}. (10) 

Thus, to evaluate fi at a given quadrature point t, only the components {Ui'}i> e s i 
need to be evaluated at t, as in Algorithm 6. This is of particular importance 
for problems of sparse structure and enables efficient multi-adaptive integration of 
time-dependent PDEs, as demonstrated below in Section 6. The sparsity pattern 
Si is automatically detected by the solver. Alternatively, the sparsity pattern may 
be specified by a (sparse) matrix. 

In Algorithm 6, the key step is the evaluation of a component Ui> at a given 
point t. For a standard mono-adaptive method, this is straightforward since all 
components use the same time steps. In particular, if the quadrature points are 
chosen to be the same as the nodal points, the value of (t) is known. For a multi- 
adaptive method, a quadrature point t for the evaluation of fi is not necessarily 
a nodal point for To evaluate Ui'(t), one thus needs to find the local interval 
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Algorithm 6 y = EvaluateRightHandSide(i, t) 

for i' E Si 

x{i')<-U v (t) 
end for 

y <- fi(x,t) 



Ii<ji such that t G Ii'j' and then evaluate Ui>(t) by interpolation on that interval. 
In Section 4 below, we discuss data structures that allow efficient storage and 
interpolation of the multi-adaptive solution. In particular, these data structures 
give 0(1) access to the value of any component in the sparsity pattern Si at 
any quadrature point t for f t . 

3.4 Multi-adaptive time step selection 

The individual and adaptive time steps kij are determined during the recursive 
construction of time slabs based on an a posteriori error estimate as discussed in 
Section 2. Thus, according to (6), each local time step kij should be chosen to 
satisfy 

k = ( tol y<«* 

where TOL is a given tolerance. 

However, the time steps can not be based directly on (11), since that leads to 
unwanted oscillations in the size of the time steps. If rij-i = max/ iJ _ 1 \Ri\ is 
small, then kij will be large, and as a result will also be large. Consequently, 
hij+i and rij + i will be small, and so on. To avoid these oscillations, we adjust 
the time step kij according to Algorithm 7, which determines the new time step 
as a weighted harmonic mean value of the previous time step and the time step 
given by (11). Alternatively, DOLFIN provides time step control based on the PID 
controllers presented in [Gustafsson et al. 1988; Soderlind 2003], including H0211 
and H211PI. However, the simple controller of Algorithm 7 performs well compared 
to the more sophisticated controllers in [Gustafsson et al. 1988; Soderlind 2003]. A 
suitable value for the weight w in Algorithm 7 is w — 5 (found empirically). 



Algorithm 7 k — Controller(fc now , fc id, k msix ) 

k <- (1 + w)k i d k ncw / (fcoid + wk new ) 
k 4- min(fc, fc max ) 



The initial time steps kn = k 2 \ = ■ ■ ■ = fcjvi = K\ are chosen equal for all 
components and are determined iteratively for the first time slab. The size K\ 
of the first time slab is first initialized to some default value, possibly based on 
the length T of the time interval, and then adjusted until the local residuals are 
sufficiently small for all components. 
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3.5 Solving the dual problem 

Stability factors may be approximated by numerically solving an auxiliary dual 
problem for (1). This dual problem is given by the following system of linear 
ordinary differential equations: 



where J(U(t),t) denotes the Jacobian of the right-hand side / of (1) at time t 
and i\) = M! (the Riesz representer of M) is initial data for the dual problem 
corresponding to the given functional Ai to be estimated. Note that we need to 
linearize around the computed solution U, since the exact solution u of (1) is not 
known. To solve this backward problem over [0,T) using the forward integrator 
Algorithm 1, we rewrite (12) as a forward problem. With w(t) — <p(T — t), we have 
w = — y>(7 — t) = J(U (T — t),T — t) T w(t), and so (12) can be written as a forward 
problem for w in the form 



For a standard mono-adaptive method, the solution on a time slab is typically 
stored as an array of values at the right end-point of the time slab, or as a list 
of arrays (possibly stored as one contiguous array) for a higher order method with 
several stages. However, a different data structure is needed to store the solution on 
a multi-adaptive time slab. Such a data structure should ideally store the solution 
with minimal overhead compared to the cost of storing only the array of degrees 
of freedom for the solution on the time slab. In addition, it should also allow for 
efficient interpolation of the solution, that is, accessing the values of the solution 
for all components at any given time within the time slab. We present below a data 
structure that allows efficient storage of the entire solution on a time slab with little 
overhead, and at the same time allows efficient interpolation with 0(1) access to 
any given value during the iterative solution of the system of discrete equations. 

4.1 Representing the solution 

The multi-adaptive solution on a time-slab can be efficiently represented using a 
data structure consisting of eight arrays as shown in Table I. For simplicity, we 
assume that all elements in a time slab are constructed for the same choice of 
method, mcG(<7) or mdG(g), for a given fixed q. 

The recursive construction of time slabs as discussed in Section 3.2 generates a 
sequence of sub slabs, each containing a list of elements (an element group). For 
each sub-slab, we store the value of the time t at the left end-point and at the right 
end-point in the two arrays sa and sb. Thus, for sub-slab number s covering the 
interval (a s ,b s ), we have 



-<j>(t) = J(U(t),t) T <p(t), te[o,T) 

<P(T) = V, 



(12) 



w(t) = f*(w(t),t) = J(U(T-t),T-t) T w(t), t e (0,T], 
tu(0) = ip. 



(13) 



4. DATA STRUCTURES 



a s = sa[s] , 
b s = sb [s] . 



(14) 
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Furthermore, for all elements in the (root) time slab, we store the degrees of freedom 
in the order they are created in the array j x (mapping a degree of freedom j to the 
value x of that degree of freedom). Thus, if each element has q degrees of freedom, 
as in the case of the multi-adaptive mcG(g) method, then the length of the array 
jx is q times the number of elements. In particular, if all components use the same 
time steps, then the length of the array jx is qN . 

For each element, we store the corresponding component index i in the array ei in 
order to be able to evaluate the correct component /j of the right-hand side / of (1) 
when iterating over all elements in the time slab to update the degrees of freedom. 
When updating the values on an clement according to (7), it is also necessary to 
know the left and right end-points of the elements. Thus, we store an array es 
that maps the number e of a given element to the number s of the corresponding 
sub-slab containing the element. As a consequence, the left end-point a e and right 
end-point b e for a given element e are given by 

a e = sa[es[e]], 
b e = sb[es[e]]. 



Array 


Type 


Description 


sa 


double 


left end-points for sub-slabs 


sb 


double 


right end-points for sub-slabs 


j x 


double 


values for degrees of freedom 


ei 


int 


component indices for elements 


es 


int 


time slabs containing elements 


ee 


int 


previous elements for elements 


ed 


int 


first dependencies for elements 


de 


int 


elements for dependencies 



Table I. Data structures for efficient representation of a multi-adaptive time slab. 



4.2 Interpolating the solution at quadrature points 

Updating the values on an element according to (7) also requires knowledge of the 
value at the left end-point, which is given as the end-time value on the previous 
element in the time slab for the same component (or the end-time value from the 
previous time slab). This information is available in the array ee, which stores for 
each element the number of the previous element (or —1 if there is no previous 
element). 

As discussed above in Section 3.3, the system of discrete equations on each time 
slab is solved by iterating over the elements in the time slab and updating the values 
on each element, either in a direct fixed-point iteration or a Newton's method. We 
must then for any given element e corresponding to some component i — ei[e] 
evaluate the right-hand side fi at each quadrature point t within the element. 
This requires the values of the solution U at t for all components contained in the 
sparsity pattern Si for component i according to Algorithm 6. As a consequence of 
Algorithm 2 for the recursive construction of time slabs, elements for components 
that use large time steps are constructed before elements for components that use 
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small time steps. Since all elements of the time slab are traversed in the same 
order during the iterative solution of the system of discrete equations, elements 
corresponding to large time steps have recently been visited and cover any element 
that corresponds to a smaller time step. The last visited element for each component 
is stored in an auxiliary array elast of size N. Thus, if i' e Si and component 
i' has recently been visited, then it is straight-forward to find the latest element 
e' = elast [z'] for component i' that covers the current element for component i and 
interpolate Uy at time t. It is also straight-forward to interpolate the values for any 
components that are present in the same element group as the current element. 

However, when updating the values on an element e corresponding to some com- 
ponent i = ei[e] depending on some other component %' e Si which uses smaller 
time steps, one must find for each quadrature point t on the element e the element 
e' for component i' containing t, which is non-trivial. The element e' can be found 
by searching through all elements for component i' in the time slab, but this quickly 
becomes inefficient. Instead, we store for each element e a list of dependencies to 
elements with smaller time steps in the two arrays ed and de. These two arrays 
store a sparse integer matrix of dependencies to elements with smaller time steps 
for all elements in the time slab. Thus, for any given element e, the number of 
dependencies to elements with smaller time steps is given by 

ed[e + 1] - ed[e], (16) 

and the elements with smaller time steps that need to be interpolated at the quadra- 
ture points for element e are given by 

{de[ed[e]],de[ed[e] + 1], . . . ,de[ed[e + 1] - 1]}. (17) 



5. PERFORMANCE 

The efficiency of multi-adaptive time-stepping compared to standard mono-adaptive 
time-stepping depends on the system being integrated, the tolerance, and the ef- 
ficiency of the implementation. For many systems, the potential speedup is large, 
but the actual speedup depends also on the overhead needed to handle the addi- 
tional complications of a multi-adaptive implementation: the recursive construction 
of time slabs and the interpolation of values within a time slab. 

To study the performance of multi-adaptive time-stepping, we consider a system 
of N components and time steps given by {kij — \Iij | : € T n } on some time slab 
T n - We define the multi- adaptive efficiency index \x by 

l/nl/'Tnax "-mm \'n\ 

where &; m i n = min/.. e 7^ kij, k max = max/,^-^ kij and \T n \ is the number of local 
intervals in the time slab T n . Thus, to obtain the multi-adaptive efficiency index, 
we divide the number of local intervals per unit time for a mono-adaptive discretiza- 
tion with the actual number of local intervals per unit time for a multi-adaptive 
discretization. This is the potential speedup when compared to a mono-adaptive 
method that is forced to use the same small time step /c m ; n for all components. 
However, the actual speedup is always smaller than yu for two reasons. The first 
is the overhead of the multi-adaptive implementation and the second is that the 
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system of discrete equations on each time slab may sometimes be more expensive 
to solve than the corresponding mono-adaptive systems (because they are typically 
larger in size). 

Consider a model problem consisting of N = Nk + Nk components, where Nk 
components vary on a slow time scale K and Nk components vary on a fast time 
scale k as in Figure 4. The potential speedup is given by the multi- adaptive effi- 
ciency index, 

= — - = — N/K — (191 

k N K + N k K/k ~~ k N K /K + N k /k ~ k ' ^ ' 

if Nk/K 3> Nk/k and K 3> k, that is the number of large elements dominates 
the number of small elements. Thus, the potential speedup can be very large for 
a system where a large part of the system varies on a large time scale and a small 
part of the system varies on a small time scale. 

If, on the other hand, K ~ k or Nk ~ Nk, then the multi-adaptive efficiency 
index may be of moderate size. As a consequence, the actual speedup may be 
small (or even "negative" ) if the overhead of the multi-adaptive implementation is 
significant. In the next section, we indicate the multi-adaptive efficiency index and 
compare this to the actual speedup for a number of benchmark problems. 



Fig. 4. A time slab with Nk = N^ = 2 and multi-adaptive efficiency index fi = 16/10 = 1.6 



6. NUMERICAL EXAMPLES AND BENCHMARK RESULTS 

In this section, we present two benchmark problems to demonstrate the efficiency 
of multi-adaptive time-stepping. Both examples are time-dependent PDEs that we 
discretize in space using the cG(l) finite clement method to obtain a system of 
ODEs, sometimes referred to as the method of lines approach. In each case, we 
lump and invert the mass matrix so as to obtain a system of the form (1). 

In the first of the two benchmark problems, the individual time steps are cho- 
sen automatically based on an a posteriori error estimate as discussed above in 
Section 3.4. For the second problem, the time steps are fixed in time and deter- 
mined according to a local CFL condition k ~ h on each element. The results were 
obtained with DOLFIN version 0.6.2. 
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6.1 A nonlinear reaction-diffusion equation 

As a first example, we solve the following nonlinear reaction-diffusion equation, 
taken from [Savcenco et al. 2005] : 



u, 



= 7u 2 (l — u) 



nx(o,n 



d n u = {) ondfix(0,T], 
«(•, 0) = uq in 17, 



(20) 



with f2 = (0, L), e — 0.01, 7 = 1000 and final time T = 1. 

The equation is discretized in space with the standard cG(l) method using a 
uniform mesh with 1000 mesh points. The initial data is chosen according to 

MX) = l + exp(A(x-l)) - (21) 

The resulting solution is a reaction front, sweeping across the domain from left to 
right, as demonstrated in Figure 5. The multi-adaptive time steps are automatically 
selected to be small in and around the reaction front and sweep the domain at the 
same velocity as the reaction front, as demonstrated in Figure 6. 
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Fig. 5. Propagation of the solution of the reaction-diffusion problem (20). 



To study the performance of the multi-adaptive solver, we compute the solution 
for a range of tolerances with L = 5 and compare the resulting error and CPU 
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x 

Fig. 6. The multi-adaptive time steps as function of space at a sequence of points in time for the 
test problem (20). 



time with a standard mono-adaptive solver that uses equal (adaptive) time steps 
for all components. To make the comparison fair, we compare the multi-adaptive 
mcG(<7) method with the mono-adaptive cG(q) method. In the benchmarks, we 
only examine q = 1. Both methods are implemented for general order q in the 
same programming language (CH — h) within a common framework (DOLFIN), but 
the mono-adaptive method takes full advantage of the fact that the time steps are 
equal for all components. In particular, the mono-adaptive solver may use much 
simpler data structures (a plain C array) to store the solution on each time slab 
and there is no overhead for interpolation of the solution. Furthermore, for the 
multi-adaptive solver, we need to supply a right-hand side function / which may 
be called to evaluate single components fi(U(t),t), while for the mono-adaptive 
solver, we may evaluate all components of / at the same time, which is usually an 
advantage (for the mono-adaptive solver). 

This is a more meaningful measure of performance compared to only measuring 
the number degrees of freedom (local steps) or comparing the CPU time against 
the same multi-adaptive solver when it is forced to use identical time steps for all 
components as in [Logg 2003a], since one must also take into account the over- 
head of the more complicated algorithms and data structures necessary for the 
implementation of multi-adaptive time-stepping. 

Note that we do not solve the dual problem to compute stability factors (or 
stability weights) which is necessary to obtain a reliable error estimate. Thus, the 
tolerance controls only the size of the error modulo the stability factor, which is 
unknown. 
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In addition, we also compare the two methods for varying size L of the domain fi, 
keeping the same initial conditions but scaling the number of mesh points according 
to the length of the domain, N = 1000L/5. As the size of the domain increases, 
we expect the relative efficiency of the multi-adaptive method to increase, since 
the number of inactive components increases relative to the number of components 
located within the reaction front. 

In Figure 7, we plot the CPU time as function of the tolerance and number of 
components (size of domain) for the mcG(l) and cG(l) methods. We also summa- 
rize the results in Table II and Table III. As expected, the speedup expressed as the 
multi-adaptive efficiency index /x, that is, the ideal speedup if the cost per degree 
of freedom were the same for the multi- and mono-adaptive methods, is large in 
all test cases, around a factor 100. The speedup in terms of the total number of 
time slabs is also large. Note that in Table II, the total number of time slabs M 
remains practically constant as the tolerance and the error are decreased. The de- 
creased tolerance instead results in finer local resolution of the reaction front, which 
is evident from the increasing multi-adaptive efficiency index. At the same time, 
the mono-adaptive method needs to decrease the time step for all components and 
so the relative efficiency of the multi-adaptive method increases as the tolerance 
decreases. See also Figure 8 for a comparison of the multi-adaptive time steps at 
two different tolerances. 

The situation is slightly different in Table III, where the tolerance is kept constant 
but the size of the domain and number of components vary. Here, the number of 
time slabs remains practically constant for both methods, but the multi-adaptive 
efficiency index increases as the size of the domain increases, since the reaction 
front then becomes more and more localized relative to the size of the domain. As 
a result, the efficiency index of the multi-adaptive method increases as the size of 
the domain is increased. 

In all test cases, the multi-adaptive method is more efficient than the standard 
mono-adaptive method also when the CPU time (wall-clock time) is chosen as a 
metric for the comparison. In the first set of test cases with varying tolerance, the 
actual speedup is about a factor 2.0 whereas in the second test case with varying 
size of the domain, the speedup increases from about a factor 2.0 to a factor 5.7 for 
the range of test cases. These are significant speedups, although far from the ideal 
speedup which is given by the multi-adaptive efficiency index. 

There are mainly two reasons that make it difficult to attain full speedup. The 
first reason is that as the size of the time slab increases, the number of iterations n 
needed to solve the system of discrete equations increases. In Table III, the number 
of iterations, including local iterations on individual elements as part of a global 
iteration on the time slab, is about a factor 1.5 larger for the multi-adaptive method. 
However, the main overhead lies in the more straightforward implementation of the 
mono-adaptive method compared to the more complicated data structures needed 
to store and interpolate the multi-adaptive solution. For constant time step and 
equal time step for all components, this overhead is roughly a factor 5 for the 
test problem, but the overhead increases to about a factor 100 when the time slab 
is locally refined. It thus remains important to further reduce the overhead of the 
implementation in order to increase the range of problems where the multi-adaptive 
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methods give a positive speedup. 




Fig. 7. CPU time as function of the error (left) and number of components N (right) for mcG(l) 
(dashed line) and cG(l) (solid line) for the test problem (20). 
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Table II. Benchmark results for mcG(l) (above) and cG(l) (below) for varying tolerance and fixed 
number of components N = 1000 for the test problem (20). The table shows the tolerance TOL 
used for the computation, the error ||e(T)||oo in the maximum norm at the final time, the time 
used to compute the solution, the number of time slabs M (with the number of rejected time slabs 
in parenthesis), the average number of iterations n on the time slab system (with the number of 
local iterations on sub-slabs in parenthesis), and the multi-adaptive efficiency index fi. 

6.2 The wave equation 

Next, we consider the wave equation, 

u tt -Au = infix(0,T], 

d n u = on9Ox(0,T], (22) 
u{-, 0) = uq in O, 
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Fig. 8. Multi-adaptive time steps at t = 0.5 for two different tolerances for the test problem (20). 
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Table III. Benchmark results for mcG(l) (above) and cG(l) (below) for fixed tolerance TOL = 
1.0T0 - 6 and varying number of components (and size of domain). (See Table II for an explanation 
of table legends.) 
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on a two-dimensional domain Q consisting of two square sub-domains of side 
length 0.5 separated by a thin wall with a narrow slit of size 0.0001 x 0.0001 at 
its center. The initial condition is chosen as a plane wave traversing the domain 
from right to left. In Figure 9, we plot the initial data together with the (fixed) 
multi-adaptive time steps. The resulting solution is shown in Figure 10. 

The geometry of the domain Q forces the discretization to be very fine close to 
the narrow slit. Further away from the slit, we let the mesh be coarse. The mesh 
was created by specifying a mesh size h with h 3> w where w is the width of the 
narrow slit. We note that for the multi-adaptive efficiency index [i defined in (19) to 
be large, the total number of elements must be large in comparison the to number 
of small elements close to the narrow slit. Furthermore, the average mesh size must 
be large compared to the mesh size close to the narrow slit. 

For a mono-adaptive method, a global CFL condition puts a limit on the size of 
the global time step, roughly given by 

k < h m i n = mmh(x), (23) 

where h = h(x) is the local mesh size. With a larger time step, an explicit method 
will be unstable or, correspondingly, direct fixed-point iteration on the system of 
discrete equations on each time slab will not converge without suitable stabilization. 

On the other hand, with a multi-adaptive method, the time step may be chosen 
to satisfy the CFL condition only locally, that is, 

k(x) < h(x), x efl, (24) 

and as a result, the number of local steps may decrease significantly (depending 
on the properties of the mesh). In this case, with k — O.lh, the speedup for the 
multi-adaptive mcG(l) method was a factor 4.2. 
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7. CONCLUSIONS 

We have presented algorithms and data structures for multi-adaptive time-stepping, 
including the recursive construction of time slabs and efficient interpolation of 
multi-adaptive solutions. The efficiency of the multi-adaptive methods was demon- 
strated for a pair of benchmark problems. The multi-adaptive methods mcG(g) 
and mdG(g) are available as components of DOLFIN, together with implementa- 
tions of the standard mono-adaptive cG(q) and dG(q) methods. The ODE solvers 
of DOLFIN are currently being integrated with other components of the FEniCS 
project, in particular the FEniCS Form Compiler (FFC) [Logg et al. 2006; Kirby 
and Logg 2006; 2007] in order to provide reliable, efficient and automatic integration 
of time dependent PDEs. 
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